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We investigate the initial-boundary value problem for linearized gravitational theory in harmonic 
coordinates. Rigorous techniques for hyperbolic systems are applied to establish well-posedness for 
various reductions of the system into a set of six wave equations. The results are used to formulate 
computational algorithms for Cauchy evolution in a 3-dimensional bounded domain. Numerical 
codes based upon these algorithms are shown to satisfy tests of robust stability for random constraint 
violating initial data and random boundary data; and shown to give excellent performance for the 
evolution of typical physical data. The results are obtained for plane boundaries as well as piecewise 
cubic spherical boundaries cut out of a Cartesian grid. 



I. INTRODUCTION 

Grid boundaries pose major difRculties in current computational efforts to simulate 3-dimensional black holes by 
conventional Cauchy evolution schemes. The initial-boundary value problem for Einstein's equations consists of the 
evolution of initial Cauchy data on a spacelike hypersurface and boundary data on a timelike hypersurface. This 
problem has only recently received mathematical attention. Friedrich and Nagy 0| have given a full solution for 
a hyperbolic formulation of the Einstein equations based upon a frame decomposition in which the connection and 
curvature enter as evolution variables. Because this formulation was chosen to handle mathematical issues rather than 
for ease of numerical implementation, it is not clear how the results translate into practical input for formulations 
on which most computational algorithms are based. The proper implementation of boundary conditions depends on 
the particular reduction of Einstein's equations into an evolution system and the choice of gauge conditions. The 
purpose of this paper is to elucidate, in a simple context, such elementary issues as: (i) which variables can be freely 
specified on the boundary; (ii) how should the remaining variables be updated in a computational scheme; (iii) how 
can the analytic results be implemented as a computational algorithm. For this purpose, we consider the evolution 
of the linearized Einstein's equations in harmonic coordinates ||^,^ and demonstrate how a robustly stable and highly 
accurate computational evolution can be based upon a proper mathematical formulation of the initial-boundary value 
problem. 

Harmonic coordinates were used to obtain the first hyperbolic formulation of Einstein's equations See Ref. ||] 
for a full account of hyperbolic formulations of general relativity. While harmonic coordinates have also been widely 
applied in carrying out analytic perturbation expansions, they have had little application in numerical relativity, 
presumably because of the restriction in gauge freedom However, there has been no ultimate verdict on the 
suitability of harmonic coordinates for computation. In particular, their generalization to include gauge source 
functions appears to offer flexibility comparable to the explicit choice of lapse and shift in conventional numerical 
approaches to general relativity. There is no question that harmonic coordinates offer greater computational efficiency 
than any other hyperbolic formulation. See Ref. j|] for a recent application to the study of singularities in a space-time 
without boundary. Here we use a reduced version of the harmonic formulation of the field equations. This allows 
us retain a symmetric hyperbolic system, and apply standard boundary algorithms, in a way that is consistent with 
the propagation of the constraints. We show, on an analytic level, that this leads to a well posed initial-boundary 
problem for linearized gravitational theory. 

Our computational results are formulated in terms of a Cartesian grid based upon background Minkowskian coordi- 
nates. Robustly stable evolution algorithms are obtained for plane boundaries aligned with the Cartesian coordinates, 
which is the standard setup for three dimensional evolution codes. Similar computational results based upon long 
evolutions with random initial and boundary data, were previously found for the linearized version of the Arnowitt- 
Deser-Misner formulation (ADM) of the field equations, where the lack of a hyperbolic formulation required 
a less systematic approach which had no obvious generalization to other boundary shapes. In this paper, we also 
attain robust stability for spherical boundaries which are cut out of the Cartesian grid in an irregular piecewise cubic 
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fashion. This success gives optimism that the methods can apphed to such problems as black hole excision and 
Cauchy-characteristic matching, where spherical boundaries enter in a natural way. 

Conventions: We use Greek letters for space-time indices and Latin letters for spatial indices, e.g. = {t, x*) 
for standard Minkowski coordinates. Linear perturbations of a curved space metric gap about the Minkowski metric 
rjajs are described by 6gaf3 — hap with similar notation for the corresponding curvature quantities, e.g. the linearized 
Riemann tensor SRap-yS and linearized Einstein tensor SGap- Indices are raised and lowered with the Minkowski 
metric, with the result that Sg"^ = —h"^. Boundaries in the background geometry are described in the form 
z = const with the spatial coordinates decomposed in the form x* = {x^,z), where the = {x,y) directions span 
the space tangent to the boundary. 



II. HARMONIC EVOLUTION WITHOUT BOUNDARIES 



The linearized Einstein tensor has the form 



where U = d^d^ 



SC^I^ = -^□7"'' - ^("r'') + ijy^'^cf^r'^, (2.1) 



^ _ (2.2) 

h = h" = -7" -7 and 

r" = -dpj"'^. (2.3) 

We analyze these equations in standard background Minkowskian coordinates x". To linearized accuracy, 

r" « .g^^V^V^x" (2.4) 

in terms of the curved space connection V^^ associated with gap and the condition F" = defines a linearized harmonic 
gauge. 

The diffeomorphisms of the curved spacetime induce equivalent metric perturbations according to the harmonic 
subclass of gauge transformations 

h°''^ ^ h^f^ + 2d^°'^'^\ (2.5) 

where the vector field f " satisfies D^" = 0. The linearized curvature tensor 

25Rf_,aui3 = d^dphai, + dad^h^_,p - dadph^^ - d^d^^hap, (2.6) 

as well as the linearized Einstein equations, are gauge invariant. 

We introduce a Cauchy foliation t in the perturbed spacetime such that it reduces to an inertial time slicing in the 
background Minkowski spacetime. The unit normal to the Cauchy hypersurfaces is given, to linearized accuracy, by 

na^-{l + lh")dat. (2.7) 

The choice of an evolution direction = an'' + with unit lapse and shift in the background Minkowski spacetime, 
defines a perturbative lapse Sa = — /i**/2 and perturbative shift (5/3* = — 



A. Standard harmonic evolution 



Harmonic evolution consists of solving the Einstein's equations Gap = 0, subject to the harmonic conditions F" = 0. 
This formulation led to the first existence and uniqueness theorems for solutions to the nonlinear Einstein equations 
by considering them as a set of 10 nonlinear wave equations [ pl| . 

In the linearized case, Einstein's equations in harmonic coordinates reduce to ten flat space wave equations so that 
their mathematical analysis is simple. The Cauchy data 7"'^(0,a:') and dt'^°'^ {0, x^) at t = determine ten unique 
solutions ^°'^{t, a;') of the wave equation 07"^ = 0, in the appropriate domain of dependence. These solutions satisfy 



2 



□r" = so that they satisfy Einstein's equations provided r"(0, x*) = and dtT°'{0, a;*) = 0, which can be arranged 
by choosing initial Cauchy data satisfying constraints. For a detailed discussion, see Ref. |l2| . 

Although this standard harmonic evolution scheme led to the first existence and uniqueness theorem for Einstein's 
equations, it is not straightforward to apply to the initial-boundary value problem. The ten wave equations for 7"^ 
require ten individual pieces of boundary data in order to determine a unique solution. Given initial data such that 
= and dtT" = at t = 0, as described above, the resulting solution satisfies the linearized Einstein equations 
only in the domain of dependence of the Cauchy data. In order that the solution of Einstein's equations extend to the 
boundary, it is necessary that F" = at the boundary. Unfortunately, there is no version of boundary data for , 
e.g. Dirichlet, Neumann or Sommerfeld data for the ten individual components, from which F" can be calculated at 
the boundary. Here we consider reduced versions of the harmonic evolution scheme in which only six wave equations 
are solved and this problem does not arise. These reduced harmonic formulations are presented below. 

B. Reduced harmonic Einstein evolution 

A linearized evolution scheme for the harmonic Einstein system, can be based upon the six wave equations 

af^ = (2.8) 

along with the four harmonic conditions 

d^r" = 0. (2.9) 

Because of the harmonic conditions, this system satisfies the spatial components 6G^^ = of the linearized Einstein's 
equations. 

As a result of the linearized Bianchi identities daSG"^ = 0, or 

dtSG'* + djSG'^ = (2.10) 
dtSG" + dj5G'^ = 0, (2.11) 

the linearized Hamiltonian constraint C := 6G** = and linearized momentum constraints := SG*^ = are also 
satisfied, throughout the domain of dependence of the Cauchy data, provided that they are satisfied at the initial 
time i = 0. This constrains the initial values of 7*" according to 

V27** + dtdjf^ ^ (2.12) 

and 

V27" - d.djf^ = 0, (2.13) 

where — S^^didj. Then, if these constraints are initially satisfied, the reduced harmonic Einstein system determines 
a solution of the linearized Einstein's equations. 

The well-posedness of the system follows directly from the wcll-posedness of the wave equations for 7*-' . The 
auxiliary variables 7"* satisfy the ordinary differential equations 

atf*+aj7''=0 (2.14) 
5*7" + 9j7*^ = 0, (2.15) 

where 7*-' enters only in the role of source terms. These differential equations do not affect the well-posedness of the 
system and have unique integrals determined by the initial values of 7"*. 

C. Reduced harmonic Ricci evolution 

The reduced harmonic Ricci system consists of the six wave equations 

□ /I'J = (2.16) 
along with the four harmonic conditions (|2.9|), which can be re-expressed in the form 
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dt(f> + djhi' + ^dth^=0, (2.17) 

+ dth'' - i + a^/i*^ = 0, (2.18) 
2 

where we have set (j) = /i**/2. Together these equations imply that the spatial components of the perturbed Ricci 
tensor vanish, SW^ = 0. In addition, the Bianchi identities imply that the remaining components satisfy 

dtC + d'C + djSR'^ = (2.19) 
dtC + dp = 0, (2.20) 

where, in terms of the Ricci tensor, C = \{6R^* + 5Rf) and C* — 5R^^. Together with the evolution equations 5R^^ ~ 0, 
the Bianchi identities imply that the Hamiltonian constraint satisfies the wave equation DC — 0. 

If the Hamiltonian and momentum constraints are satisfied at the initial time, then dtC also vanishes at the initial 
time so that the uniqueness of the solution of the wave equation ensures the propagation of the Hamiltonian constraint. 
In turn, Eq. ( p.l9D then ensures that the momentum constraint is propagated. Thus the reduced harmonic Ricci system 



of six wave equations and four harmonic equations leads to a solution of the linearized Einstein's equations for initial 
Cauchy data satisfying the constraints. 

The harmonic Ricci system takes symmetric hyperbolic form when the wave equations are recast in first differential 
order form. Thus the system is well posed. The formulation of the harmonic Ricci system as a symmetric hyperbolic 
system and the description of its characteristics is given in Appendix Although well-posedness of the analytic 
problem does not the guarantee the stability of a numerical implementation it can simplify its attainment. 



D. Other reduced harmonic systems 

The harmonic Einstein and Ricci systems are special cases of a one parameter class of reduced harmonic systems 
for the variable 

which satisfying the wave equations E^^ = —^CIk^^ = and the harmonic conditions = 0, where E"^ = SG"^ — 
^rj^^SG. This system is symmetric hyperbolic when the auxiliary system F" = is symmetric hyperbolic. This can 
be analyzed by setting k*^ = in the auxiliary system, which then takes the form 

= a..--^^aVS (2.21) 

and implies 

A 



{dt - J^x^f^d^y = 0- (2-22) 



The auxiliary system is symmetric hyperbolic when Eq. ( 2.22 ) is a wave equation whose wave speed 



A 



(3A-2) 



(2.23) 



is positive. This is satisfied for A < and A > 2/3. In the range 2/3 < A < 1, the wave speed is faster than the speed 
of light. Only for the case A = 1, the harmonic Ricci system, is the wave speed of the auxiliary system equal to the 
speed of hght. 

The auxiliary system for the reduced harmonic Einstein case has a well posed initial-boundary value problem but 
represents a borderline case. This could adversely affect the development of a stable code based upon a nonlinear 
version of the reduced harmonic Einstein system. 
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III. THE INITIAL-BOUNDARY VALUE PROBLEM 



Consider the initial-boundary problem consisting of evolving Cauchy data prescribed in a set S at time t — 
lying in the half-space z > with data prescribed on a set T on the boundary z = Q. Harmonic evolution takes its 
simplest form when Einstein equations are expressed as a second differential order system. However, in order to apply 
standard methods it is necessary to recast the problem as a first order symmetric hyperbolic system. Then the theory 
determines conditions on the boundary data for a well posed problem in the future domain of dependence . Here 
is the maximal set of points whose past directed inextendible characteristic curves all intersect the union of S 
and T before leaving I?+ . 

Appendix ^ describes the symmetric hyperbolic description of boundary conditions for the 3-dimensional wave 
equation. The basic ideas and their application to the initial-boundary value problem are simplest to explain in terms 
of the 1-dimensional wave equation, as follows. See Appendix |^ for the analogous treatment of the 3-dimensional 



case. Our presentation is based upon the formulation of maximally dissipative boundary conditions |13|, the approach 
used by Friedrich and Nagy in the nonlinear case. An alternative description of boundary conditions for symmetric 
hyperbolic systems is given in Ref . and for linearized gravity in Ref . [ 15| . 

A. The 1-dimensional wave equation 

The one-dimensional wave equation h^tt — h^zz can be recast as the first order system of evolution equations 

ii^k (3.1) 
k = (3.2) 
/ = (3.3) 

where we have introduced the auxiliary variables k and /. Given initial data (/i, k, /) at time t = subject to 
the constraint, / — h z, these equations determine a solution h of the wave equation. Note that the constraint is 
propagated by the evolution equations. 



In coordinates = {x^,x^) = {t,z) the system (3.1)-( p^ ) has the symmetric hyperbolic form 

A^'df^u ^Bu + F, (3.4) 
where the solution u consists of the column matrix 

/ h\ 

w = I fc I (3.5) 

and 

/10 0\ /000\ /010\ 

A" = 010 , ^1= 0-1, B= 000 . (3.6) 

yooiy \o-ioy yoooy 

In our case the source ter m F — but otherwise it plays no essential role in the analysis of the system. 
The contraction of Eq. (3.4) with the transpose give the flux equation 



^ uA'^df.u ^'u{'B + B)u. (3.7) 

This can be used to provide an estimate on the norm J ^uA^udz for establishing a well posed problem, in the half-space 
z > with boundary at z = 0, provided the flux arising with the normal component of A^ satisfies the inequality 

'^uA^u < 0. (3.8) 

This inequality determines the allowed boundary data for a well posed initial-boundary value problem. 

As expected from knowledge of the characteristics of the wave equation, the normal matrix A^ has eigenvalues 
0,±1, with the corresponding eigenvectors 
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^^0 = , «±i = Tl . (3.9) 



These eigenvectors are associated with the variables 

ho^h, h±i^ fTk^h,,Th,t. (3.10) 
Re-expressing the solution vector in terms of the eigenvectors, 





(3.11) 



the inequality (3.8 implies that homogeneous boundary data must take the form 

= h+i- Hh^i, (3.12) 

where the parameter H satisfies < 1. The component ho, corresponding to the kernel of A^, propagates directly 
up the boundary and cannot be prescribed as boundary data. 
Non-homogeneous boundary data can be given in the form 

q{t) = h+i - Hh^i = /i,^(l -H)- h^t{l + H), (3.13) 

where q is an arbitrary function representing the free boundary data at z = 0. Eq. ( 3.13| ) shows how the scalar 
wave equation accepts a continuous range of boundary conditions. The well-known cases of Sommerfeld, Dirichlet or 
Neumann boundary data are recovered by setting H ~ 0, -|-1, or, —1, respectively. Note that there are consistency 
conditions at the edge {t = 0, z — 0). For instance, Dirichlet data corresponds to specifying k = h^t on the boundary 
and this must be consistent with the initial data for k. 



B. The reduced harmonic Einstein system 



The harmonic Einstein system consists of the six wave equations (2.8) and the four harmonic conditions ( p.Sj ). 
Since the wave equations for 7*-' are independent of the auxiliary variables 7*", the well-posedness of the initial- 
boundary value problem for 7'^ follows immediately. Furthermore, since the harmonic conditions propagate the 
auxiliary variables up the boundary by ordinary differential equations, the harmonic Einstein system has a well posed 
initial-boundary value problem. A unique solution in the appropriate domain of dependence is determined by the 
initial Cauchy data 7'^ and 7^ at t = 0, the initial data 7*" at the edge t = z = and the boundary data 7'-' at z = 
given in any of the forms described in Appendix ^ (e.g. Dirichlet, Neumann, Sommerfeld). 

A solution of the linearized harmonic Einstein system satisfies SG^^ =0. As a result the Bia nchi identit ies im ply 



dtC^ = and dtC = —djC^ so that the constraints are satisfied provided the constraint Eq's. ( 2.1^ ) and ( 2.15 ) are 
satisfied at t — 0. 

The free boundary data for this system consists of six functions. However, as shown in Ref. [|l|, the vacuum Bianchi 
identities satisfied by the Weyl tensor imply that only two independent pieces of Weyl data can be freely specified 
at the boundary. We give the corresponding analysis for the linearized Einstein system in Appendix This makes 
it clear that only two of the six pieces of metric boundary data are gauge invariant. This is in accord with the four 
degrees of gauge freedom consisting of the choice of linearized lapse Sa = — /i*'/2 (one free function) and linearized 
shift 5/3* = — /i*' (three free functions). 

A linearized evolution requires a unique lapse and shift, whose values can be specified explicitly as space-time 
functions or specified implicitly in terms of initial and boundary data subject to dynamical equations. In the case of 
a harmonic gauge, in order to assess whether explicit space-time specification of the lapse and shift is advantageous 
for the purpose of numerical evolution, it is instructive to see how it affects the initial-boundary value problem for the 
reduced harmonic Einstein system. In the linearized harmonic formulation without boundary, a gauge transformation 
( ^.5|) to a shift-free gauge 7** = is always possible within the domain of dependence of the initial Cauchy data . 
In the presence of a boundary, consider a gauge transformation with ^" = (0,^*), so that 

=7"-f-a'f . (3.14) 
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For harmonic evolution of constrained initial data, both and 7** satisfy the wave equation. At t = 0, choose Cauchy 
data for ^* satisfying 

= 7" + 9*e (3.15) 

and 

= da*' + dtd*C = ~dj-f'' ~ V^f . (3.16) 
On the boundary z — 0, require that ^* satisfy 

= 7*^+a*r. (3.17) 

Then 7** = 7*' + 9*^* has vanishing Cauchy data at t = 0, vanishing Dirichlet boundary data at z = and satisfies 
the wave equation, so that 7** — 0. Thus, for evolution of constrained data with a boundary, a shift-free harmonic 
gauge is possible. 

However, in this gauge, boundary data for all 6 components of 7*'' can no longer be freely specified since the 
harmonic condition implies 

9.7""' + 5^7^^ = 0. (3.18) 

This relates Neumann data for to Dirichlet data for and would complicate any shift-free numerical evolution 
scheme. As an example, one could freely specify Dirichlet boun dary data for the 3 components 7^^^ and (i) obtain 
Neumann boundary data dz^^^ from the A-components of Eq. ( |3.18 ), (ii) evo lve 7 ^"^ in terms of initial Cauchy data 



and (iii) obtain Neumann boundary data 9x7^^ from the z-componcnt of Eq. ( |3.18 ). Note the nonlocality of step (ii) 



which would have to be carried out "on the fly" during a numerical evolution. 

The requirement that the shift vanish reduces the free boundary data to three components. It is also possible to 
eliminate an additional free piece of data by choosing a unit lapse, i.e. setting /i** = 0. Suppose the shift has been 
set to zero, so that the harmonic condition implies dth** = —dth\. Consider a gauge transformation /i** — /i** — 29t^*, 
where ^" — (f*,C*) satisfies dtC — <9*^* so that the shift remains zero. For harmonic evolution of constrained data, 
both ^* and /i" satisfy the wave equation. At t = 0, choose Cauchy data for ^* satisfying 



and 



= /i" - 2dti' (3-19) 



= dth}* - 2diC = -dth\ - 2V^e (3.20) 



(where we assume the Cauchy data is given on a non-compact set so that there are no global obstructions to a solution 
^*). On the boundary z = 0, require ^* satisfy 

= /i"-2c)t^*. (3.21) 

Then /i** — because it is a solution of the wave equation with vanishing Cauchy data and Dirichlet boundary data. 
(Alternatively, the lapse can be gauged to unity by a transformation satisfying F* = D^* = const, so that the harmonic 



source function F* still drops out of Eq. (2.1) for the Einstein tensor.) 

A unit lapse and zero shift implies that d^l = so that 7- cannot be freely specified at the boundary. Coupled with 
our previous results, imposition of a unit lapse and zero-shift reduces the free boundary data to the two trace-free 
transverse components — (l/2)^^^7p, in accord with the two degrees of gauge-free radiative freedom associated 



with the Weyl tensor. (In addition, the initial Cauchy data must satisfy Eq. (3.18), dtjl — and V^7- = 0). A 
similar result arises in the study of the unit-lapse zero-shift initial-boundary value problem for the linearized ADM 
equations However, in the case of harmonic evolution, it is clear that explicit specification of the lapse and shift 
leads to a more complicated initial-boundary value problem. It is more natural to retain the freedom of specifying 6 
pieces of boundary data, which then determine the lapse and shift implicitly during the course of the evolution. 

The reduction of the free boundary data can be accomplished by other gauge conditions on the boundary, which are 
not directly based upon the lapse and shift. An example, which plays a central role in the Friedrich-Nagy formulation, 
is the specification of the mean curvature of the boundary. To linearized accuracy, the unit outward normal to the 
boundary at z = is 

A^a « (1 + i/i..)V„z. (3.22) 
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The associated mean extrinsic curvature (g"^ — N" N^)'V isNa is given to linear order by 

X = dtht. - dAhf + \{d,hi ~ d^hu). (3.23) 

Although the extrinsic curvature of a planar boundary vanishes in the background Minkowski space, the linear 
perturbation of the background induces a non-vanishing linearized extrinsic curvature tensor. 
Under a gauge transformation induced by the the mean curvature transforms according to 

X = x + (5,'-5^5^)r. (3.24) 

A gauge deformation of the boundary in the embedding space-time makes it possible to obtain any mean curvature 
by solving a wave equation intrinsic to the boundary. In this respect, the mean extrinsic curvature of the boundary is 
pure "boundary gauge" and can be specified to eliminate one degree of gauge freedom. When the harmonic conditions 
are satisfied, the mean curvature of the boundary reduces to x = ^9^/1^^. 

This discussion shows that there are various ways that the six free pieces of boundary data can be restricted by 
gauge conditions. Such restrictions can be important for an analytic understanding of the initial-boundary problem 
but their usefulness for numerical simulation is a separate issue, especially in applications where the boundary does 



not align with the numerical grid as discussed in Sec. IV 



C. The reduced harmonic Ricci system 



The underlying eq uation s of t he re duced harmonic Ricci system are the six wave equations (2. If) and the four 
harmonic conditions ( 2.17 ) and ( 2.1g ). The symmetric hyperbolic formulation of this system and the analysis of its 
characteristics is given in Appendix |A|. The variables consist of 

{h'\T'\X'\Y'\Z'\h'\<j)) , 



where T*-' = dth^^ , X^^ = d^h^^ , Y"^^ — dyh^^ and Z'-* = d^K^K In a well posed initial-boundary value problem, there 
are seven free functions that may be specified at the boundary. For example, in the analogue of the Dirichlct case, 
the free boundary data consists of T'-' and (j). 



In the context of the second order differential form given by Eq's. ( 2.16 ), with the harmonic conditions ( 2.17 ) and 
(2.18), the boundary data remains the same, e.g. Dirichlet data dth"^^ and 0. This determines th e evo lution of K^ ^ 
by the wave equation, which then provides source terms for the symmetric hyperbolic subsystem ( ^.17 ) and ( 2.18 ). 
This subsystem implies that = 0, so that the evolution of (j) is also governed by t he w ave equation, with initial 



Cauchy data <j) and dt4> where dt4> is provided by the initial Cauchy data /i** via Eq. (2.18)). The evolution of ft.** is 



then obtained by integration of Eq. ( p.l8| ). 

Unlike the reduced harmonic Einstein system where the constraints propagate up the boundary by ordinary dif- 
ferential equations, the initial-boundary value problem for the reduced harmonic Ricci system d oes n ot necessarily 
satisfy Einstein's equations even if the constraints are initially satisfied. The Bianchi identities ( ^.19 ) imply ( 2.20| ) 
so that both C and dtC would initially vanish but, since C satisfies the wave equation, C would vanish throughout 
the evolution domain only if it vanished on the boundary. In that case, Eq. (2.19) would imply that the momentum 
constraints were also satisfied throughout the evolution domain. 

Thus evolution of constrained initial data for the harmonic Ricci system yields a solution of the Einstein equations 
if and only if the Hamiltonian constraint is satisfied on the boundary. This is equivalent to requiring that = on 
the boundary. If the evolution equations are satisfied, we can express this in the form 



3 * 



0, 



(3.25) 



where H — h\ and H"^^ = h^^ — i(5*-'7i^. This allows formulation of the following well posed initial-boundary value 
problem for a solution satisfying the constraints. 

We prescribe initial Cauchy data that satisfies the constraints for i/*-', H, <t> and /i'*, and free b oundary data for 
H^^ and (j). The system of wave equations then determines H''^ . This allows integration of Eq. ( 3.25 ) on the boundary 
to obtain Dirichlet boundary values for determination of as a solution of the wave equation. The rem ainin g fields (f> 
and /i** are then determined as a symmetric hyperbolic subsystem. Note that the boundary constraint ( 3.25 ) reduces 
the free boundary data from seven independent functions (for unconstrained solutions) to six, in agreement with the 
free boundary data for solutions of the reduced harmonic Einstein system. 
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IV. NUMERICAL IMPLEMENTATION 



Numerical error is an essential new factor in the computational implementation of the preceding analytic results. 
The initial Cauchy data cannot be expected to obey the constraints exactly. In particular, machine roundoff error 
always produces an essentially random component to the data which, in a linear system, evolves independently of 
the intended physical simulation. It is of practical importance that a numerical evolution handle such random data 
without producing exponential growth (and without an inordinate amount of numerical damping). We designate 
as robustly stable an evolution code for the linearized initial-boundary problem which does not exhibit exponential 
growth for random (constraint violating) initial data and random boundary data. This is the criterion previously 
used to establish robust stability for ADM evolution with specified lapse and shift 

We test for robust stability using the 3-stage methodology proposed for evolution-boundary codes in Ref. The 
tests check that the norm of the Hamiltonian constraint C does not exhibit exponential growth under the following 
conditions. 



Stage I: Evolution on a 3-torus with random initial Cauchy data. 



• Stage II: Evolution on a 2-torus with plane boundaries, i.e. x [— L, L], with random initial Cauchy data and 
random boundary data. 

• Stage III: Evolution with a cubic boundary with random initial Cauchy data and random boundary data. 

Stage I tests robust stability of the evolution code in the absence of a boundary. Stage II tests robust stability of 
the boundary-evolution code for a smooth boundary (topology T^). Stage III tests robust stability of the boundary- 
evolution code for a cubic boundary with faces, edges and corners, as standard practice for computational grids based 
upon Cartesian coordinates. 



We have established robust stability for the evolution-boundary codes, described in See's. IV A and IVB, the 
reduced harmonic Einstein and Ricci systems. The tests are performed according to the procedures outlined in Ref. 
1^ for an evolution of 2000 crossing times (a time of 4000L, where 2L is the linear size of the computational domain) 
on a uniform 48'^ spatial grid with a time step At = Aa;/4 (which is slightly less than half the Courant-Friedrichs-Lewy 
limit and typical of values used in numerical relativity) . 

For purposes such as singularity excision or Cauchy-characteristic matching, there are computational strategies 
based upon spherical boundaries. The extension of a robustly stable ADM boundary-evolution algorithm with given 
lapse and shift from a cubic grid boundary to a spherical boundary is problematic fl^ . Robustly stable evolution- 
boundary algorithms for the reduced harmonic Ricci and Einstein systems with spherical boundaries are presented in 



Sec. [VC 



Nimrerical evolution of the fields (or 7'^) and ft,** (or 7**) is implemented on a uniform spatial grid 
{^[i]:y[j]7 Z[K]) = {I J Ax, KAx), with time levels i'^l = N At. Thus a field component / is represented 
by its values flf\ j,^ = /(il^l i^[i]iy[j]^Z[K])- In order to obtain compact finite difference stencils for impos- 
ing boundary conditions, the fields /i** (or 7'*) are represented on staggered grids at staggered time-levels, where 

^^-1-1/2] j-/4-fW+l/2l \ 

J[I+l/2,J+l/2,K+l/2] ' ' , X[i+i/2],y[J+l/2], Z[K+l/2])- 

The extensive tests of stability reported here were performed using a leapfrog evolution algorithm in order not 
to bias the tests by introducing excessive dissipation. We expect that the tests would also be satisfied by more 
dissipative algorithms. This was borne out by a limited number of evolutions using an iterative Crank-Nicholson 
algorithm (implemented as in Ref. 0). 



A. Robustly stable algorithms for the reduced harmonic Einstein system 

Stage I 

Evolution of 7'-' is done by a standard three-level leapfrog scheme, with the wave equation in second-differential- 
order form, i.e. 

VIN+1]_^ vim vlN-l] 
^[I,J,K] '^^[I,J,K]^^[IJ,K] __tj[N] ,^VIN] ,^^m] (A^\ 

(^At)'^ ~ '^.^2;[7,J,if] l,yylI,J,K] l,zzlI,J,K] ' V^--^; 

Second spatial derivatives are computed via 

r _ /[/+!. J.A'] - f[I,J.K] + f[I-l.J,K] , s 

J,xx[I.J,K] - • l4.^j 
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On a g rid st aggered in space and at staggered time levels, evolution of 7'* is carried out by the finite difference version 
of Eq. (I2II ), i.e. 



it[N+l/2] _ it\N-l/2] 

^[I+l/2,J+l/2,K+l/2\ ^[I+\/2,J+l/2,K+l/2\ _ 

At 

ixlN] iy[N] izlN] 

~'^,x[I+l/2.,J+l/2,K+l/2] ~ '^,y[I+l/2,J+l/2,K+l/2] ~ 'T,z[/+l/2, J+l/2,if+l/2] ' 

Here first spatial derivatives are evaluated at the center of the integer grid cells, i.e. 



(4.3) 



1 



flI+l,J,K] — flI,J,K] /[/+!,, 7+1, if] — flIJ+l,K] 



J,x[I+l/2,J+l/2,K+l/2\ — ^ XZ ^ 7^ ^ > 



Ax Ax 



^ /[/+!, J.jf+l] ~ f[I,J,K+l\ _|_ f[I+l,J+lJ<+l] — f[I,J+lJ<+l] 



Ax Ax 

Since 7** is represented on the integer grid which is staggered in space and ti me w ith respect to the grid representation 
of 7'*, the finite difference equation used for updating 7" is similar to Eq. (O) used to update 7**. 

Stage II 

The hierarchy of ten linearized equations for the reduced harmonic Einstein system makes the initial-boundary value 
problem particularly simple to implement as a finite difference algorithm. Let the plane boundary be given by the 
-fCo-th grid point, with the Stage I evolution algorithm applied to all points with K < Kq. Update of 7[*^^_]^] requires 

7^^^j as boundary data. The same boundary information allows update of Y[Ko-i/2]^ which, in turn, allows update of 
"/[Ko-i] ^ith no additional boundary data. However, it is interesting to note that specification of free boundary data 
for 7*^jj] or 7*^q_|_]^/2] '^o^ld not affect stability simply because no evolution equation uses that data. 

Stage III 

Let the evolution domain be the cube — L < x* < L with the grid structured so that the boundary lies on (non- 
staggered) grid points, i.e., x|-^j = ±L, etc. The boundary data consist of 7*^ on the faces, edges and corners of the 
cube. The field 7'* can then be updated at all staggered grid points inside the cube, including those neighboring the 
boundary. For instance, update of 7p^2 3/2 3/2] involves use of at the points [| ± i, | ± i, | ± i], all of which are 
on or inside the boundary. Similarly, evolution of 7*' can be carried out at all interior grid points without further 
boundary data. 

Robust stability of the evolution-boundary algorithm is demonstrated by the graph of the Hamiltonian constraint 
in Fig. 0. The linear growth results from momentum constraint violation in the initial data. 



B. Robustly stable algorithms for the reduced harmonic Ricci system 



Stage I 



Evolution of h^^ is carried out identically as the evolution of 7*.^ in the Einstein system (see Eq. (4T)). The fields 
h^^ and are represented on the int eger g rid w hile / t'* is represented on a half-integer grid staggered in space and in 
time. Thus the evolution equations ( 2.17) and (2.18) for ip and ft.** have finite difference form 



dm 



At 



lN+1/2] 



At 



MN+l/2] 



MN-1/2] 



"[7+1/2, .7+1/2, AT -H/2] '"[/+1/2,. 7+1/2, X+1/2] 

At 



a.ft" - -Sjkd'h^' 



[N] 



= 0, 



(4.5) 
(4.6) 



[7+1/2,, 7+1/2, 7f+l/2] 



where the spatial derivative terms are computed according to Eq. (4.4) 

Stage II 



Unconstrained plane boundary. As shown in Sec. IIIC, seven free functions can be prescribed as unconstrained 



boundary data for h^^ and cf) in the reduced harmonic Ricci system. Again let the boundary be defined by the Kq-IYi 
grid point, with K < for interior points. Then the boundary data h^^Xo] ^'^'^ '^[^0] allows the evolution algorithm 
to be appfied to update h^^ and /i** at all interior points, e.g. to update h^^i^^_^-^ and ft[Xo_i/2]' which in turn allows 
update of (j) at all interior points. 
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Constrained plane boundary. Conservation of the constraints in the reduced harmonic Ricci system requires that 
the Hamiltonian constraint be enforced at the boundary. In order to obtain a finite difference approximation to a 
solution of Einstein's equations, the unconstrained evolution-boundary algorithm must be modified to enforce the 
Hamilt onian constraint on the boundary. On a z =const boundary, we accomplish this, in accord with the discussion 
in Sec. [lie, by prescribing freely the function (f) and the five components of t he traceless symmetric tensor dtH'^-' ■ 
The missing ingredient, dtH is updated at the boundary according to E g. (3.25 ). 

In the finite difference algorithm, in order to be able to apply Eq. ( 3.25| ) via a centered three-point stencil, we 
introduce a guard point at z = zj/f^^]^], where the boundary data dtH^-' is provided. (See Fig. |l|.) Assuming that the 

fields /i*^ , ft.'' and (p are known for t < t^^\ K < Kq and that the boundary data H^^^^_^_^ is known for t < t^^\ we 

use the following evolution-boundary algorithm to compute these fields at tl^+^l; 

• (i) We update the fields h^^ at time-level tl^+^l at all grid points within the numerical domain of dependence of 
the data known at t^^l, i.e., at all points which require no boundary data. 



• (ii) At the guard point Kg -I- 1 we assign Dirichlet boundary data to the five independent components of the 
symmetric traceless (ft^'"') |^|^(]^' by prescribing boundary values for dtH'^^ , dtH^^ , dtH^^ , dtH^y , dtH^^ and 
setting = -dtH'''' - dtRyy . 



• (iii) At the guard point /Cq + 1 we update the fields H^^^'^^-^ 



[7V+1/2] 
'[^^0 + 1] ■ 



• (iv) At the boundary point Kq we update H^^^^^ 
difference form 



i][N+l] 



using the field-equation □i?*-' = 0, written in the finite- 



rrij[N+l] _ r, TTij[N] , TTij[N-l] 

Ai2 



(4.7) 



• (v) At the boundary point Kq we compute the boundary values H^^^^^ using the finite-difference version of 



Eq. ( |3.25D , i.e 



2 ^[A-o] 



1] 9 TriN] 



At2 



[N-l] 
[Ko] 



Am 



0. 



(4.8) 



In Eq's. (4/7) - (4^) all derivatives are computed using centered three-point expressions. 



• (vi) From knowledge of iJ'^'^^""^' and h]^^^'^'^ we construct /i*^^^"'^'. 



(vii) We assign boundary data for 



and update ^I^+^l and /i»*[^+i/2] according to Eq's. (gj) 



(4.6) 



Stage III 



Unconstrained cubic boundary. The unconstrained cubic boundary is essentially the same as the unconstrained 
plane boundary, i.e., the functions K^^ and are provided at the boundary. Then the evolution equations are applied 
to update ft'-' , ft** and cj) at all interior grid points. Recall that ft*' is represented on a staggered grid with some 
inter ior p oints located a half grid-step away from the boundary. Nevertheless, ft'* can be updated at such points using 
Eq. (Ul). 

Constrained cubic boundary. The algorithm for enforcing the Hamiltonian constraint on a cubic boundary is an 
extension of the algorithm for a plane boundary. The boundary functions dtH''-' are provided at a set of guard points 
that are Ax outside the boundary of the cube, such that at grid-points on the boundary of the cube we can use the 
field equations □i/*-' = as well as the boundary constraint Eq. ( [3.25 ) to update and H at the boundary points. 
Boundary data for is provided at grid points on the boundary of the cube, in a similar fashion to the constrained 
plane boundary. 

Robust stability of the evolution-boundary algorithm for the constrained cubic boundary is demonstrated by the 
graph of the Hamiltonian constraint in Fig. ^. In addition to the robust stability test results, in Fig. |^ we have also 
included results from two physical runs, based on the plane- wave solution Eq's. ( 1.11 ) - ( 1.14 ). The longer physical 
run (up to t/L « 1000) was performed with a grid-size of Ax — 1/80, while the shorter run (up to t/L « 400) was 
performed with a grid-size of Ax = 1/120. Comparison of the physical runs with different gridsizes indicates that the 
long-term polynomial growth in the Hamiltonian constraint violation can be controlled by grid-resolution. 
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C. Spherical boundaries 



The implementation of both the reduced harmonic Einstein and Ricci systems display robust stability in Stage I, II 
and III tests. We now extend these results to a a spherical boundary cut out of a Cartesian grid, with the evolution 
domain defined by the interior of a sphere of radius R. In order to update fields at these interior grid points, values 
are needed at a set of "guard" points consisting of grid points on or outside the boundary. Some of these field values 
constitute free boundary data. We extend our test of robust stability to a fourth stage which checks that the ioo norm 
of the Hamiltonian constraint C does not exhibit exponential growth under the following condition. 

• Stage IV: Evolution with a spherical boundary with random initial Cauchy data and random free boundary 
data at all guard points. 

As before, we perform the tests for an evolution of 2000 crossing times (4000i?, where 2R is the diameter of the 
computational domain) on a uniform 48'^ spatial grid with a time step slightly less than half the Courant-Friedrichs- 
Lewy limit. We have established Stage IV robust stability for the following implementations of reduced harmonic 
Einstein and Ricci revolution. 



1. The reduced harmonic Einstein system 

Let the evolution domain be defined by the interior of a sphere of radius R. The update of 7** at an interior 
point [I,J,K] requires the values of 7** at the points [/± ^j^/i ^^-^^ il' some of which might be "staggered- 
boundary" points outside the evolution domain. These staggered-boundary points are inside the spherical shell 
R < ^/x^ + y^+z^ < R + ^ 6rQ, if Svq > ^/3Ax. At these points we update 7** by the same finite-difference equation 



(4.3) as for the interior points. This defines a set of non-staggered boundar y po ints for the field 7*^ determined by two 
conditions: (i) that the 7*^ can be updated at all interior points using Eq. (4T) and (ii) that 7** can be updated at all 



interior and staggered-boundary points. Both conditions are satisfied if we define the set of non-staggered boundary 
points by i? < W^^^/j + vfj^ + ^^^j < R + Sro, where Stq > \/3 Ax. At this set of points we obtain values for 7*-' via 



Eq. (4.1), with 9t7*^ provided in a "guard-shell" R + Stq < ^Jx?^^^%p^^f^z^ < R + Sro + Svi, where Svi > Ax. The 
radius R of the spherical boundary of the reduced harmonic Einstein system is related to the linear size 2L of the 
computational grid by 

R + 5ro + Sri<L. (4.9) 
This ensures that all guard points fall inside the domain [—L, +L]'^. 



2. The reduced harmonic Ricci system 



Unconstrained spherical boundary. The Stage IV unconstrained evolution-boundary algorithm for the reduced 
harmonic Ricci system is similar to that of the reduced harmonic Einstein system. We use two spherical shells, of 
thickness (5ro and (5ri, with 6rQ > ^/SAx and Sri > Ax. The algorithm is the following: 

• (i) We update all fields h^^ at time- level ^I^+^l at all non-staggered grid points within the sphere of radius 
R + 6ro. 

• (ii) We provide boundary data {dthY''^^^'^^^^ at all non-staggered grid points within the spherical shell R+dro < 
^/x^ + y^ + z'^ < R + Sro + dri. We update = h'^^^^+At {dthf^^^^^^^ within the same spherical shell. 



• (iii) We provide boundary data 01^+^1 within the spherical shell R < x^ + + < R + SrQ. 

• (iv) We update the field inside the sphere of radius R ac cordi ng to Eq. ( ^.17 ), and update the fields 
^it[jv+i/2] ij^gi(jg ^i^g sphere of radius R + (5ro/2 according to Eq. ( 2.18 ). 

Constrained spherical boundary. The evolution-boundary algorithm for the reduced harmonic Ricci system with 
constrained spherical boundary is an extension of the algorithm with unconstrained boundary. In addition to the 
spherical shells defined in the case of the unconstrained boundary, we define an additional set of guard points R + 
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'in 



y'^j-^ + zj^^j < i? + (5ro + Sri + Sr2, where the boundary data dfH'^-' is provided. The quantity Sr2 



Sra + Sri < 

is defined by two conditions: (i) the fields H^^ can be updated according to Eq. (4.7) at all non-staggered grid points 
R + Sro < ^J^Ji] + y[j] + ^[K] < + '^'^o + '^'^i, and (ii) the Hamiltonian constraint Eq. ( 3.25 ) can be enforced at the 

same set of grid points. Both of these conditions are satisfied if Sr2 > Ax. The evolution-boundary algorithm for 
the reduced harmonic Ricci system with a constrained spherical boundary is the following: 

• (i) We update all fields h^^ at time-level tl^+^l at all non-staggered grid points within the sphere or radius 
R + Sro. 



(ii) We update the fields iJ'Jt^+i] at the set of boundary points R + Sro < 

via Eq. Then we update the field i/I^+il at these points via Eq. (^ 

we construct /i^^I^+i] at the same set of grid points. 

(iii) We provide boundary data '^"""^^^ 

R + Sro + Sri + Sr^, then we update iJ^I^+i] = + M {ptH 

(iv) We assign boundary data for at the set of boundary points R < 

(v) We update ^I^+^l and /i**[^+i/2] according to Eq's. - ( p^ . 



Sro + Sri 
From knowledge of and 



^[/] + y[j] 



at the set of guard points R - 



V] 



Sro + Sri < 
at the same set of grid points 

Sro. 



'[I] ^ y[j] 



yf.n + ^K] < R 



The radius R of the constrained spherical boundary for the reduced harmonic Ricci system and the linear size 2 L 
of the computational grid are related by 



R + Sro + Sri + Sr2 < L. 



(4.10) 



This ensures that all boundary points and all guard points fall inside the domain [—L, -l-L]"^. 

The graphs of the Hamiltonian constraint in Fig. ^ illustrate robust stability for a spherical boundary for the reduced 
harmonic Einstein system and for the reduced harmonic Ricci system with constrained boundary data. Comparison 
of Figs. H - ^ shows that there is no significant difference between the Stage HI and Stage IV performances in terms 
of numerical stability. 



D. Convergence Tests 

In order to calibrate the performance of the algorithms we carried out convergence tests based upon analytic 
solutions constructed from a superpotential <I>"'^^'' symmetric in (a,/3) and (/i,i^) and antisymmetric in [a,//], [/?, i^], 
such that □$"'3'^'' = 0. As a result of these symmetry properties, the tensor 7^"^ = dadp^°'^^'^ is symmetric and 
satisfies the linearized harmonic Einstein equations □7'^'^ — and d^'y^^ = 0. 

In our first testbed we choose as a superposition of two solutions, 

^A.M^'^ = A^"" + B'"', fi^i^, (4.11) 
with the remaining independent components of ^^'^^'^ set to zero. The solution A'^'^ is defined by 

A" = A'* = A** = 0, (4.12) 

^,^8«-^-sin[..(^-Ja.-+^-)/^]^ (4.13) 



and i?^'' is defined by 



8 6^*sinra;_B (^-x*)] 

^ L_|_V 11^ = B" = 0. (4.14) 



Here A^^ is a plane wave propagating with frequency u>a along the diagonal of the(a;, y) plane, so that a wave 
crest leaving x = —L travels a distance before arriving at a; = +L. Since the topology of Stages I and II imply 
periodicity in the x direction, we set 



UJA 



= V2tt/L. 
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Similarly, the frequency Wb of the frinctions S** is set to 

oJb = tt/L. 

In Stage III we use the same choices, while in the Stage IV tests we set 

UJA — V^tt/R, lob = tt/R- 

The amplitudes a^^ , 6** were chosen to be 

a^'' = 1.1 • 10-^ a^^ = 1.3•10-^ a^'^ = 1.2 • 10-^ 
6^* = 1.4•10-^ 6*'* = 1.0 • 10-^ 6^* = 1.5-10-^ 

Convergence runs used the plane wave solution. In Stages I - III we used the grid sizes 

Ax _ 1 1 1 1 
2L ~ 80' Too' 120' 160 ■ 

while in Stage IV we used 

Aa; _ 1 1 1 1 
2^ ^ 80' lOO' 120' 160 ' 

with the additional gridsize Ax/{2R) = 1/200 in the case of the the reduced harmonic Ricci system with constrained 
spherical boundary. 

The time-step was set to At = Ax /A. In the Stage IV test of the reduced harmonic Einstein system the widths of 
the boundary shells were chosen to be Stq = 1.8 Ax and 6ri = Ax. The same parameters were used when testing the 
algorithm for the reduced harmonic Ricci system with unconstrained spherical boundary. The evolution-boundary 
algorithm for the reduced harmonic Ricci system with constrained spherical boundary was tested using the parameters 
Sro = 1.8 Ax, dri = 2.5Ax, and dr2 = 1.5 Ax. 

The code was used to evolve the solutions from t = to t = L {t = R in. Stage IV), at which time convergence 
was tested by measuring the ioo and the £2 norms of ^7** for the Einstein system and of □(^ for the Ricci system, 
which test convergence of the Hamiltonian constraint. The norms were evaluated in the entire evolution domain. In 
addition, we also checked convergence of the metric components to their analytic values. 

In addition to plane wave tests, we tested the qualitative performance in Stage IV using an offset spherical wave 
based upon the superpotential (with shifted origin) 

^,,..^ fit + r)-nt-f) ^ (4.15) 

where f = ^/x'^ + {y + a)^ + and 

The parameters A, a, and w are set to 

A = 2-w^ ■ 10"^ a = OMR, w = 0.2R. 



-( 



uy 
wJ 



(4.16) 



1. The reduced harmonic Einstein system 

Evolution requires Cauchy data at t = and boundary data at the guard points. The Cauchy data 
{7*-' , 9(7*-' , 7**, 7**}j^g was provided by giving -y^ilo] ^ ^ijl-i] ^ ^ttlo] g^^^^ ^«*[i/2] gj[ interior and guard points. In 

addition, we provided boundary data at each time-step by giving (9*7'^)'^^^^^' at all guard points. The metric and 
Hamiltonian constraint were found to be 2nd order convergent for Stages I - IV. In particular, in Stage IV, the norm 
of 07" vanished to 0{A^-^^). 
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2. The reduced harmonic Ricci system 



In the case of the Ricci system the Cauchy data { ft,'-' , dt ft'-' , ft'* , 0} was provided by giving ft'-' [''I , ft*^ [ ^1 , (Z)!"! and 
/j«*[i/2] q]\ interior and boundary points. In addition, when the Hamiltonian constraint was numerically imposed 
at the boundary, we also provided by giving iJ^^I"!. 

We first tested the code without numerically imposing the Hamiltonian constraint. In this case we provided 

boundary data at each time-step by giving (atft^^)"^+'^'l and at all guard points. 

Next we tested the code with the Hamiltonian constraint numerically imposed at the boundary. Thus we first 
prescribed the traceless ( dtH^ ^) t^^-'-^^l ^^d (/>[^1 at all guard points, then computed H at each time-step via the 
boundary constraint Eq. ( [3. 251 ). 

In all cases we found the numerically evolved metric functions converge to their analytic values to O(A^). In Stage 
I, vanished to roundoff accuracy, while in Stages II - III it vanished to second order accuracy. In particular, for 
the Stage III algorithm with constrained boundary, dip converged to zero as 0(A^-^^). 

In Stage IV, with constrained boundary, we found that the £2 norm of Ocj) vanishes to first order accuracy. However, 
the £00 norm decreases linearly with grid size only for Ax/{2R) = 1/80,1/100 and 1/120 but fails to show further 
decrease for Ax/{2R) — 1/160 and 1/200. This anomalous behavior of the £00 norm stems from the random way in 
which guard points are required at different sites near the boundary. This introduces an unavoidable nonsmoothness 
to the second order error in the metric components, which in turn leads to 0(1) error in the second spatial derivatives 
occurring in □</> or in the Hamiltonian constraint. Unlike the Einstein system in which the constraints propagate 
tangent to the boundary, this error in the Ricci system propagates along the light cone into the interior. However, 
since its origin is a thin boundary shell whose width is 0{Ax), the £2 norm of remains convergent to first order. We 
expect that the convergence of the Hamiltonian constraint for a spherical boundary would be improved by matching 
the interior solution on the Cartesian grid to an exterior solution on a spherical grid aligned with the boundary, as is 
standard practice in treating irregular shaped boundaries. 



3. Simulation of an outgoing wave using a constrained spherical boundary 



We also tested the code's ability to evolve an outgoing spherical wave traveling off center with respect to a spherical 
boundary of radius R. Fig. ^ illustrates a simulation performed using the Stage IV algorithm for the reduced harmonic 
Ricci system, with the Hamiltonian constraint numerically enforced at the boundary. The metric fields were evolved 
from t = to t = 1.5R, using a grid of Ax/{2R) = 1/120. After the analytic wave has propagated out of the 
computational domain, the remnant error is two orders of magnitude smaller than the initial signal. This shows that 
artificial refiection off the boundary is well controlled even in the computationally challenging case of a piecewise cubic 
spherical boundary. 
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APPENDIX A: HYPERBOLIC FORMULATION OF THE LINEARIZED HARMONIC RICCI SYSTEM 



In order to study the evolution of the system consisting of Eq's. (2.16), ( ^.17 ), and (2.18) in the half-space z > with 
boundary at z = 0, we employ the auxiliary variables T*-' = dth^^ ^X^^ = dxh^-' ,y^-' = dyh^^ ,Z'^^ = dzh^-',4> 
In terms of the variables 



ift" 



the system takes the form 
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^2 2 2 

dthy' = -dy(i> + ^Y'"'' ~ ^Yyy + ^ y"- - x-^y - z^^ 
dth^* = -d^cj) + ^z'''= + ^zyy z^^ - x'^'- - y^" 

dt<\, = -d^h""^ - dyhy' - d,h'' --T''''--Tyy-- T'-\ 

X y z 2 2 2 

Next we define the 34-dimensional vector u by 

u = "^(/i^^, h'^y, /^^^ hyy, h^, /i^", t^^, . . . , t^-, x^^, . . . , x'', r^^^, . . . , z^^, ... ,z 



The system of equations ( |AlD - (A9) then has the form 
where ^* = I34X34 is the identity matrix, 



A'^d^u = Bu 



and 



where 
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The matrix has the eigenvalue +1, with muhiphcity 7 and eigenvectors 
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the eigenvalue —1, with multiphcity 7 and eigenvectors 
and the kernel of the matrix has dimension 20, with a basis 



K'KhyK 



In the eigen-basis defined by A'' the vector u defined in Eq. (AlO) takes the form 



with 



uo 

u+ 



^XX rj^ZZ 



J^J22 YXX 

Z'-\6 + h'-*). 



Non-homogeneous boundary data can be given in terms of a free column vector field q in the form 

q = — 

where H can be any 7x7 matrix satisfying 



aa- 



(A18) 



(A19) 
(A20) 
(A21) 



(A22) 



(A23) 



The three simplest matrices that satisfy the condition ( A23 ) are —1^x7, Oyx? and lyx? . The first of these corresponds 
to specifying Neumann data for h^^ and Dirichlet data for /i^*. Using the zero matrix as a candidate for H corresponds 
to giving Sommerfeld data for h^^ and specifying the quantity </> — /i^*. Last, picking H to be the identity matrix 
corresponds to giving Dirichlet data for dth''-' as well as for </>. Note that the evolution system (Al) - (A9) accepts 
a much richer clas s of boundary conditions than the three we just mentioned. One simply needs to pick a matrix H 
that satisfies Eq. (A23) and the choice of H defines the seven free functions that are to be specified at the boundary. 



APPENDIX B: WEYL DATA ON BOUNDARY 

The curvature tensor, which provides gauge invariant fields, decomposes into the Ricci curvature, which vanishes 
if the evolution and constraint equations are satisfied, and the Weyl curvature Cap-yS- In order to analyze the 
boundary freedom, it is convenient make the following choice of a complete, independent set of 10 linearized Weyl 
tensor components Kap^s = SCap-yS- Cab = KtABt - ^SabS'^^ KtcDt, Kt^At, KtAzu Kt^xy, Kt^AB, and Dab = 
KtABz — ^^abS^cdz- linearized vacuum Bianchi identities d^^Kp^^^^, — and d^KapjS = to show that 

the Weyl data which can be freely specified on the boundary can be reduced to the 2 independent components C'ab ■ 

First, the identity d^KtzzS = implies (after using the trace-free property of the Weyl tensor) 



dtK, 



't^t At 



ab 



^ 



which determines the boundary behavior of Ct'^At in terms of the remaining 9 Weyl components. 
Next, note that the identity 



dtKtABC + dBKtACt + dcKtAtB = 0. 



implies 



dtKt^Ac = -d^CAc + +::dcKt''Du 



or, taking a i-derivative and using Eq. (Bl), that 



dfRt'^Ac = -dtd^CAc + Idcd'^Kt^AD 



(Bl) 



(B2) 



(B3) 



(B4) 
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This gives a propagation equation intrinsic to the boundary which determines the time dependence of Kt^AC in terms 
of the boundary data for Cab- (Note that Kt^AC propagates up the boundary with velocity in one mode and in a 
cone with velocity l/\/2 in the other mode.) 
Next, the identity 

d'KtAzS - -d'KtAzt + d^'KtAzB = (B5) 

determines the time dependence of CtAzt', and the identity 

dtKtzAB + dBKtztA + dAKtzBt = (B6) 

determines the time dependence of Ktzxy This reduces the free Weyl data on the boundary to the 4 independent 
components Cab and Dab- 

However, the specification of Dab, in addition to Cabi would lead to an inconsistent boundary value problem. 
This can be seen from the identity 

d^KtABS = -OtKtABt + d.KtABz + d'^KtABC = 0, (B7) 

which determines Neumann data for Dab in terms of Dirichlet data for Cab and other known quantities. Similarly, 
the identity 

dtKtABz - d.KtABt + dBKtAzt = (B8) 

determines Neumann data for Cab- Thus, since the components of the Weyl tensor satisfy the wave equation, the 
specification of both Cab and Dab as free Dirichlet boundary data leads to an inconsistent boundary value problem. 

The determination of boundary boundary values for Dab from boundary data for Cab is a global problem which 
first requires solving the wave equation to determine Cab from its boundary and initial data. Then the time derivative 
of the trace-free part of Eq. (B8) yields 

O'^^Dab - ObO^Dac + ^^ABd^'d^'DcD - dtd^CAB - (B9) 



which propagates Dab up the boundary in terms of initial data. Defining C = q^q^Cab and D = q^q^ Dab, with 
q-^Oa = dx + idy, this reduces to 

dfD ~ ^OaO^D - dtd.C = 0, (BIO) 

which has propagation velocity 1 / V2. (Note that this is but one of the variations consistent with the maximally 
dissipative condition used by Friedrich and Nagy Q]. In the case of unit lapse and vanishing shift, assigning boundary 
data for C is equivalent to assigning data for the trace-free part of the intrinsic 2-metric of the boundary foliation, 
consistent with results found in Ref. |0].) 
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[Ko-2] [Kq-I] [Kq] [Kq+I] 

FIG. 1. Boundary stencil for the Stage 2 evolution-boundary algorithm of the Ricci system with constrained boundary. 
Circles stand for interior grid points which are updated by the evolution algorithm. The required boundary data for H^^ is 
provided at the guard point [Kq + 1] (triangle), while the boundary data for is provided at the boundary point [Ko\ (square). 
The Hamiltonian constraint is enforced at the boundary point [Kq\. 
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FIG. 2. A log-log plot of the l^o norm of the Hamiltonian constraint as a function of time for a Stage 3 test of the 
evolution-boundary algorithm of the Einstein system and of the Ricci system with constrained boundary. The upper two 
curves correspond to stability tests (random data), while the lower two curves indicate performance tests (physical data). 
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FIG. 3. A log-log plot of the iao norm of the Hamiltonian constraint as a function of time for a Stage 4 test of the 
evolution-boundary algorithm of the Einstein system and the Ricci system with constrained boundary. 
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FIG. 4. Stage IV evolution of an outgoing solution using the reduced harmonic Ricci system with a constrained spherical 
boundary. The plots show the metric field /i"' at a; = 0, for t/ R — 0.1 (top), t/R = 1.0 (middle) and t/R = 1.5 (bottom). In 
the bottom plot, the field has decayed by two orders of magnitude. 
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